#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Feb 17 11:44:57 2024

@author: jcfq2
"""

import os

jwst_dir='/Users/jcfq2/data/observations/jwst'

os.chdir(jwst_dir)

import matplotlib.pyplot as plt
import matplotlib.colors as colours
from astropy.visualization import make_lupton_rgb
from astropy.table import Table
import ch4_fiddlesticks as ch4
import pandas as pd
import h3ppy
import glob
import spiceypy as spice
from astropy.io import fits
import numpy as np
from importlib import reload
import JWSTSolarSystemPointing as jssp
from scipy.interpolate import interp1d

reload(jssp)
#import jwst_uranus as jwstu

# In[0]

kernel_dir = '/Users/jcfq2/data/observations/jwst/kernels/'
jssp.load_kernels(kdir=kernel_dir)

# Set up a h3ppy object too, always useful
h3p_model = h3ppy.h3p()

# In[0] set up new colormaps



colors2 = plt.cm.Greens_r(np.linspace(0, 1, 220))
colors1 = plt.cm.cubehelix_r(np.linspace(0.0, 0.4, 32))

# combine them and build a new colormap
colors = np.vstack((colors1, colors2))

for i in range(32):
    colors[31-i,0]=colors[32,0]*(32-(i+1))/32
    colors[31-i,1]=colors[32,1]*(32-(i+1))/32
    colors[31-i,2]=colors[32,2]*(32-(i+1))/32
# for i in range(32): print(colors[32,1]*(32-(i+1))/32)

mygreen = colours.LinearSegmentedColormap.from_list('my_green', colors)


colors2 = plt.cm.Reds_r(np.linspace(0, 1, 220))
colors1 = plt.cm.cubehelix_r(np.linspace(0.0, 0.4, 32))

# combine them and build a new colormap
colors = np.vstack((colors1, colors2))

for i in range(32):
    colors[31-i,0]=colors[32,0]*(32-(i+1))/32
    colors[31-i,1]=colors[32,1]*(32-(i+1))/32
    colors[31-i,2]=colors[32,2]*(32-(i+1))/32
# for i in range(32): print(colors[32,1]*(32-(i+1))/32)

myred = colours.LinearSegmentedColormap.from_list('my_red', colors)


colors2 = plt.cm.Blues_r(np.linspace(0, 1, 220))
colors1 = plt.cm.cubehelix_r(np.linspace(0.0, 0.4, 32))

# combine them and build a new colormap
colors = np.vstack((colors1, colors2))

for i in range(32):
    colors[31-i,0]=colors[32,0]*(32-(i+1))/32
    colors[31-i,1]=colors[32,1]*(32-(i+1))/32
    colors[31-i,2]=colors[32,2]*(32-(i+1))/32
# for i in range(32): print(colors[32,1]*(32-(i+1))/32)

myblue = colours.LinearSegmentedColormap.from_list('my_blue', colors)



colors2 = plt.cm.Purples_r(np.linspace(0, 1, 220))
colors1 = plt.cm.cubehelix_r(np.linspace(0.0, 0.4, 32))

# combine them and build a new colormap
colors = np.vstack((colors1, colors2))

for i in range(32):
    colors[31-i,0]=colors[32,0]*(32-(i+1))/32
    colors[31-i,1]=colors[32,1]*(32-(i+1))/32
    colors[31-i,2]=colors[32,2]*(32-(i+1))/32
# for i in range(32): print(colors[32,1]*(32-(i+1))/32)

mypurple = colours.LinearSegmentedColormap.from_list('my_purple', colors)


img = np.arange(350**2).reshape((350, 350))

plt.imshow(img,cmap=myblue)


# In[9] as above but with four movie panels
from JWSTSolarSystemPointing import get_pixel_polygons as jssp_gpp

from pypolyclip import clip_multi, clip_single
 

xscale=100
yscale=100

#450 for 30 degrees of wrap around either size
# naxis is weird, if it is too small it skips anything outside the range, but can be any size bigger than the range of values clipped

naxis = (10*xscale,10*yscale)

naxis_mapcount=np.zeros((int(10*xscale+1),int(10*yscale+1)))
naxis_h3p=np.zeros((int(10*xscale+1),int(10*yscale+1)))
naxis_methane=np.zeros((int(10*xscale+1),int(10*yscale+1)))
naxis_heat=np.zeros((int(10*xscale+1),int(10*yscale+1)))
naxis_reflect=np.zeros((int(10*xscale+1),int(10*yscale+1)))
naxis_F323=np.zeros((int(10*xscale+1),int(10*yscale+1)))
naxis_lat=np.zeros((int(10*xscale+1),int(10*yscale+1)))
naxis_lon=np.zeros((int(10*xscale+1),int(10*yscale+1)))
naxis_ra=np.zeros((int(10*xscale+1),int(10*yscale+1)))
naxis_dec=np.zeros((int(10*xscale+1),int(10*yscale+1)))
files = sorted(glob.glob(
    "/Users/jcfq2/data/observations/jwst/5308_dither_separated/*.fits"))

# files=files[10]

lens=len(files)
# lens=5
z0=0
zz=0
x_offset=5
y_offset=-3

savedir="/Users/jcfq2/OneDrive - Northumbria University - Production Azure AD/python/jwst/saturn/movie2/"


# minorlontick=[0,15,30,60,75,90,105,120,150,165,180,195,210,240,255,270,285,300,330,345]
# majorlontick=[45,135,225,315]


majorlontick=[0,90,180,270]
minorlontick=[45,135,225,315]


for i in range(lens):
    
    
    
    # file = 'data/jw03665032001_02101_g395h-f290lp_s3d.fits'
    geo = jssp.JWSTSolarSystemPointing(files[i])
    wave = geo.get_wavelength()
    cube = geo.full_fov()
    spec = geo.convert(wave, geo.im[:, 25, 25])
    cube = geo.full_fov()
    cube_1 = geo.full_fov(corner=1)
    cube_2 = geo.full_fov(corner=2)
    cube_3 = geo.full_fov(corner=3)
    cube_4 = geo.full_fov(corner=4)
    

    # build a wave corrected F323N filter
    removeH3p = 'True'
    
    if i == 0:
        # correction for 3953
        wavemin_F323 = 3.148
        wavemax_F323 = 3.35
        whw_F323 = np.argwhere((wave > wavemin_F323) & (wave < wavemax_F323)).flatten()

        F323_values=np.loadtxt('F323N_filter.txt',delimiter=' ')

        f_F323 = interp1d(F323_values[:,0], F323_values[:,1], kind='linear')  # You can change 'linear' to other interpolation methods like 'cubic'
        filter_F323_sub = f_F323(wave[whw_F323])

        filter_F323 = np.zeros_like(wave)
        filter_F323[whw_F323]=filter_F323_sub
        filter_F323_h3p = filter_F323*1.0
        filter_F323_without_h3p = filter_F323*1.0
        
        
        h3p = h3ppy.h3p()
        h3p.set(T=500, N=6e14, wave=wave, R=2000)
        h3p_model=h3p.model()
        h3p_model=h3p_model/np.max(h3p_model)
        
        filter_F323_h3p[h3p_model < 0.001]=0
        filter_F323_without_h3p[h3p_model > 0.001]=0
    
        # filter_F323 = filter_F323_h3p

    sat_lat = cube[0, :, :]
    sat_lon = cube[1, :, :]
    ra = (cube[13, :, :]-geo.ra_target)* -3600.0         *xscale +(x_offset*xscale)
    ra_1 = (cube_1[13, :, :]-geo.ra_target)* -3600.0     *xscale +(x_offset*xscale)
    ra_2 = (cube_2[13, :, :]-geo.ra_target)* -3600.0     *xscale +(x_offset*xscale)
    ra_3 = (cube_3[13, :, :]-geo.ra_target)* -3600.0     *xscale +(x_offset*xscale)
    ra_4 = (cube_4[13, :, :]-geo.ra_target)* -3600.0     *xscale +(x_offset*xscale)
   
    dec = (cube[14, :, :]-geo.dec_target)* 3600.0       *yscale +(y_offset*yscale)
    dec_1 = (cube_1[14, :, :]-geo.dec_target)* 3600.0   *yscale +(y_offset*yscale)
    dec_2 = (cube_2[14, :, :]-geo.dec_target)* 3600.0   *yscale +(y_offset*yscale)
    dec_3 = (cube_3[14, :, :]-geo.dec_target)* 3600.0   *yscale +(y_offset*yscale)
    dec_4 = (cube_4[14, :, :]-geo.dec_target)* 3600.0   *yscale +(y_offset*yscale)


    sy=3
    ey=-2
    geoim=geo.im[:, :,sy:ey]
    sat_lat = sat_lat[:,sy:ey]
    sat_lon = sat_lon[:,sy:ey]
    ra = ra[:,sy:ey]
    dec=dec[:,sy:ey]
    ra_1 = ra_1[:,sy:ey]
    dec_1=dec_1[:,sy:ey]
    ra_2 = ra_2[:,sy:ey]
    dec_2=dec_2[:,sy:ey]
    ra_3 = ra_3[:,sy:ey]
    dec_3=dec_3[:,sy:ey]
    ra_4 = ra_4[:,sy:ey]
    dec_4=dec_4[:,sy:ey]
    wavemin = 3.0
    wavemax = 3.2
    
    # dlambda=0.00015
    #whw = np.argwhere((wave > 3.3529+dlambda) & (wave < 3.3535+dlambda)).flatten()
    whw = np.argwhere((wave > wavemin) & (wave < wavemax)).flatten()
    img_reflect = np.nanmedian(geoim[whw, :, :], axis=0)
    
    img_methane = geoim[2775, :, :]
    
    
    aaa=np.nanmean(geoim[672:677,:,:],axis=0)
    bbb=np.nanmean(geoim[672+7:677+7,:,:],axis=0)
    # ccc=np.nanmean(nam[:,:,672+22:677+22]/nac[:,:,672+22:677+22],axis=2)
    # ddd=np.nanmedian(geoim[whw_331,:,:],axis=0)

    img_methane = geoim[2775, :, :]
    img_methane = bbb-aaa

    
    wavemin = 3.9529
    wavemax = 3.9535
    
    whw = np.argwhere((wave > wavemin) & (wave < wavemax)).flatten()
    img_h3p = np.nanmedian(geoim[whw, :, :], axis=0)
    
    ab=geoim[2239,:,:]+geoim[1568,:,:]+geoim[1642,:,:]+geoim[1010,:,:]+geoim[1011,:,:]+geoim[1018,:,:]+geoim[1019,:,:]
    cd=(geoim[2241,:,:]*0.8*0.5+geoim[2236,:,:]*1.19*0.5)+((geoim[1571,:,:]+geoim[1565,:,:])*0.5)+(0.4*geoim[1640,:,:]/5.21)+(geoim[1004,:,:]+geoim[1005,:,:]*2)

    img_h3p = ab-cd
    
    img_heat = geoim[3487,:,:]

    #ra, dec = geo.get_delta_lat_lon_arcsec()
    
    # im_h3p[im_h3p > np.nanmedian(im_h3p)*3] = np.nanmedian(im_h3p)*3
    # im_h3p[im_h3p<0] = 0
    
    i0=0

    for xx in range(ab[:,0].size):
        for yy in range(ab[0,:].size):
        
            
            
            pixel_h3p=img_h3p[xx,yy]
            pixel_methane=img_methane[xx,yy]
            pixel_heat=img_heat[xx,yy]
            pixel_reflect=img_reflect[xx,yy]
            
            pixel_F323=np.nansum(geoim[:,xx,yy]*filter_F323)
            # print(xx,yy,pixel_int)
            pixel_lat = sat_lat[xx,yy]
            pixel_lon = sat_lon[xx,yy]
            pixel_ra = ra[xx,yy]
            pixel_dec = dec[xx,yy]
            
            
            if ~np.isnan(pixel_h3p):
                
                pxx = np.array([ra_1[xx,yy],ra_2[xx,yy],ra_3[xx,yy],ra_4[xx,yy]])
                pyy = np.array([dec_1[xx,yy],dec_2[xx,yy],dec_3[xx,yy],dec_4[xx,yy]])
                pii_h3p=pixel_h3p
                pii_methane=pixel_methane
                pii_heat=pixel_heat
                pii_reflect=pixel_reflect
                pii_F323=pixel_F323
                pii_lat=pixel_lat
                pii_lon=pixel_lon
                pii_ra=pixel_ra
                pii_dec=pixel_dec


                if i0 == 0:

                    px=pxx
                    py=pyy
                    pi_h3p=pii_h3p
                    pi_methane=pii_methane
                    pi_heat=pii_heat
                    pi_reflect=pii_reflect
                    pi_F323=pii_F323
                    pi_lat=pii_lat
                    pi_lon=pii_lon
                    pi_ra=pii_ra
                    pi_dec=pii_dec
                else: #i0 is set to one on first iteration, if not set, then make px and py

                    px= np.vstack([px, pxx])
                    py= np.vstack([py, pyy])
                    pi_h3p= np.hstack([pi_h3p, pii_h3p])
                    pi_methane= np.hstack([pi_methane, pii_methane])
                    pi_heat= np.hstack([pi_heat, pii_heat])
                    pi_reflect= np.hstack([pi_reflect, pii_reflect])
                    pi_F323= np.hstack([pi_F323, pii_F323])
                    pi_lat= np.hstack([pi_lat, pii_lat])
                    pi_lon= np.hstack([pi_lon, pii_lon])
                    pi_ra= np.hstack([pi_ra, pii_ra])
                    pi_dec= np.hstack([pi_dec, pii_dec])
                    
                i0=i0+1
                    # print(xx,yy,pxx,pyy,pii)

                # print(i0,pxx,pyy)

    xc, yc, area, slices = clip_multi(px, py, naxis)
 
#     # compute the total area by summing over all the pixels for each polygon

    A0 = np.asarray([np.sum(area[s]) for s in slices], dtype=float)
 
#     # make arrays to capture the values

    m_counts=np.zeros_like(area)
    m_counts_noadds=np.zeros_like(area)
    m_h3p=np.zeros_like(area)
    m_methane=np.zeros_like(area)
    m_heat=np.zeros_like(area)
    m_reflect=np.zeros_like(area)
    m_F323=np.zeros_like(area)
    m_lat=np.zeros_like(area)
    m_lon=np.zeros_like(area)
    m_ra=np.zeros_like(area)
    m_dec=np.zeros_like(area)

#     # itrator for slices/pixels (where slices are bunched by pixels)

    ixel=0

    for s in slices: 

        m_counts[s] = area[s]#/np.sum(area[s])
        m_counts_noadds[s]= 1

        # pypolyclip gives a value that is the proportion of the lat long pixel (which might be right?), here, I instead use the proportion of the pixel, so all latlong subpixels sum to a count of 1 - might be the wrong way to do this, but biases towards the smallest (most face on) pixels.  Needs tetsing!
        m_h3p[s] = pi_h3p[ixel]*area[s]#(area[s]/np.sum(area[s]))
        m_methane[s] = pi_methane[ixel]*area[s]#(area[s]/np.sum(area[s]))
        m_heat[s] = pi_heat[ixel]*area[s]#(area[s]/np.sum(area[s]))
        m_reflect[s] = pi_reflect[ixel]*area[s]#(area[s]/np.sum(area[s]))
        m_F323[s] = pi_F323[ixel]*area[s]#(area[s]/np.sum(area[s]))
        m_ra[s] = pi_ra[ixel]*area[s]#(area[s]/np.sum(area[s]))
        m_dec[s] = pi_dec[ixel]*area[s]#(area[s]/np.sum(area[s]))
        m_lat[s] = pi_lat[ixel]#(area[s]/np.sum(area[s]))
        # m_lon[s] = pi_lon[ixel]*area[s]#(area[s]/np.sum(area[s]))
        m_lon[s] = pi_lon[ixel]#(area[s]/np.sum(area[s]))
        ixel=ixel+1
 
    # fill individual map positions with the above

    if z0 < 4:
        if z0 == 0:
            start_time = geo.obs_start
            cml = (geo.lon_sun+360)%360
            geofirst=geo
       
    # naxis_mapcount[yc,xc]=naxis_mapcount[yc,xc]+m_counts
    # naxis_h3p[yc,xc]=naxis_h3p[yc,xc]+m_int
        naxis_mapcount[yc,xc]=m_counts+naxis_mapcount[yc,xc]
        naxis_h3p[yc,xc]=m_h3p+naxis_h3p[yc,xc]
        naxis_methane[yc,xc]=m_methane+naxis_methane[yc,xc]
        naxis_heat[yc,xc]=m_heat+naxis_heat[yc,xc]
        naxis_reflect[yc,xc]=m_reflect+naxis_reflect[yc,xc]
        naxis_F323[yc,xc]=m_F323+naxis_F323[yc,xc]
        naxis_lat[yc,xc]=m_lat
        naxis_lon[yc,xc]=m_lon
        naxis_ra[yc,xc]=m_ra+naxis_ra[yc,xc]
        naxis_dec[yc,xc]=m_dec+naxis_dec[yc,xc]
        
        z0=z0+1
        
        
        
        
    else:
        z0=1
        zz=zz+1
        image_g = np.nan_to_num(naxis_h3p/naxis_mapcount)
        image_g=np.nan_to_num(image_g)**0.5
        image_g=image_g/135*3
        image_g[image_g != 0]= image_g[image_g != 0]-0.2
        image_g[image_g < 0]=0
        
        image_r = np.nan_to_num(naxis_heat/naxis_mapcount)  # /np.max(im_methane)
        # print(np.nanmax(image_r))
        # image_r=image_r/np.max(image_r)
        # image_r=image_r/140000*1.1
        image_r=image_r/140000*3
        image_F323 = np.nan_to_num(naxis_F323/naxis_mapcount)  # /np.max(im_reflect)
        image_reflect = np.nan_to_num(naxis_reflect/naxis_mapcount)  # /np.max(im_reflect)
        image_b = np.nan_to_num(naxis_methane/naxis_mapcount)  # /np.max(im_reflect)
        # print(np.nanmax(image_b))
        # image_b=image_b/np.max(image_b)
        # image_b=image_b/90000*5
        image_b=image_b/90000*150
        image = make_lupton_rgb(image_r*5, image_g*5, image_b*10)


        print(zz,np.nanmax(naxis_lon))

        # naxis_lon=naxis_lon2
        # image_lat=np.nan_to_num(naxis_lat)
        # image_lon=np.nan_to_num(naxis_lon)
        
        image_ra=np.nan_to_num(naxis_ra/naxis_mapcount)
        image_dec=np.nan_to_num(naxis_dec/naxis_mapcount)
        
        image_ra2=  ((image_ra-(x_offset*xscale))/xscale/-3600.0)+geofirst.ra_target
        image_dec2= ((image_dec-(y_offset*yscale) )/yscale/3600.0)+geofirst.dec_target

        image_ra2[image_ra == 0]=0
        image_dec2[image_dec == 0]=0
        
                
        
        image_lat=np.zeros_like(image_ra)*np.nan
        image_lon=np.zeros_like(image_ra)*np.nan
        
        
        for xx in range(image_ra[:,0].size):
            for yy in range(image_ra[0,:].size):
            
                
                if image_ra2[xx,yy] != 0:
                    geopixel=geofirst.pixel_params(image_ra2[xx,yy],image_dec2[xx,yy])
                    image_lat[xx,yy]=geopixel['lat']
                    image_lon[xx,yy]=geopixel['lon']
       
        # if image_ra2[xx,yy] != 0:
        #     rara=((yy-(x_offset*xscale))/xscale/-3600.0)+geofirst.ra_target
        #     decdec=((xx-(y_offset*yscale) )/yscale/3600.0)+geofirst.dec_target
        #     geopixel=geofirst.pixel_params(rara,decdec)
        #     image_lat[xx,yy]=geopixel['lat']
        #     image_lon[xx,yy]=geopixel['lon']

        if np.abs(geofirst.lon_sub) < 90:
            image_lon2=image_lon+0.
            image_lon2[image_lon2>180]=image_lon2[image_lon2>180]-360
            image_lon2[np.nonzero(image_lon2)]=image_lon2[np.nonzero(image_lon2)]+180
            image_lon=image_lon2+0.


        # image_lat2=image_lat+0
        # image_lon2=image_lon+0
    # for xxx in range(180): aaa[:,xxx]=aaa[:,xxx]/np.nanmax(aaa[:,xxx])

                
        fig = plt.figure(figsize=(16,9))
        
        scaling = ((700*3/16)/(460*2/9))
        # xwidth*xnum*16/yheight*ynum*9
        
        mytext=start_time[0:10]+' '+start_time[11:19]
        
        a0 = plt.subplot(position=[0,1.0/scaling,1,0.164])
        a0.set_axis_off()
        a0.text(0.5, 0.7, mytext,c='w',stretch='ultra-expanded',family='Arial Black', ha="center", va="center", weight='bold',size=24)

        a1 = plt.subplot(position=[0,0.5/scaling,1/3,0.5])
        # a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
        a1.set_axis_off()        
        
        a1.imshow(image,aspect=1,origin='lower')
        # plt.title(str(i))
        # plt.plot([0,180],[550,550],linestyle='--',color='grey')
        a1.set_xlim((1*xscale,8*xscale))
        a1.set_ylim((0.9*yscale,5.5*yscale))
        
        a1.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
        a1.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
        a1.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)
        
        a2 = plt.subplot(position=[1/3,0.5/scaling,1/3,0.5])
        # a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
        a2.set_axis_off()        
        
        a2.imshow(image_g,aspect=1,origin='lower',vmin=0,cmap=mygreen)
        # plt.title(str(i))
        # plt.plot([0,180],[550,550],linestyle='--',color='grey')
        a2.set_xlim((1*xscale,8*xscale))
        a2.set_ylim((0.9*yscale,5.5*yscale))
        
        a2.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
        a2.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
        a2.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)

        
        a3 = plt.subplot(position=[2/3,0.5/scaling,1/3,0.5])
        # a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
        a3.set_axis_off()        
        
        a3.imshow(image_r,aspect=1,origin='lower',vmin=0,cmap=myred)
        # plt.title(str(i))
        # plt.plot([0,180],[550,550],linestyle='--',color='grey')
        a3.set_xlim((1*xscale,8*xscale))
        a3.set_ylim((0.9*yscale,5.5*yscale))
        
        a3.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
        a3.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
        a3.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)

        
        a4 = plt.subplot(position=[0,0.0,1/3,0.5])
        # a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
        a4.set_axis_off()        
        
        a4.imshow(image_reflect,aspect=1,origin='lower',vmin=0,cmap='Greys_r')
        # plt.title(str(i))
        # plt.plot([0,180],[550,550],linestyle='--',color='grey')
        a4.set_xlim((1*xscale,8*xscale))
        a4.set_ylim((0.9*yscale,5.5*yscale))

        a4.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
        a4.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
        a4.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)

        a5 = plt.subplot(position=[1/3,0.0,1/3,0.5])
        # a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
        a5.set_axis_off()        
        
        a5.imshow(image_b,aspect=1,origin='lower',vmin=0,cmap=myblue)
        # plt.title(str(i))
        # plt.plot([0,180],[550,550],linestyle='--',color='grey')
        a5.set_xlim((1*xscale,8*xscale))
        a5.set_ylim((0.9*yscale,5.5*yscale))


        a5.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
        a5.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
        a5.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)

        a6 = plt.subplot(position=[2/3,0.0,1/3,0.5])
        # a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
        a6.set_axis_off()        
        
        a6.imshow(image_F323,aspect=1,origin='lower',vmin=0,cmap='copper')
        # plt.title(str(i))
        # plt.plot([0,180],[550,550],linestyle='--',color='grey')
        a6.set_xlim((1*xscale,8*xscale))
        a6.set_ylim((0.9*yscale,5.5*yscale))

        a6.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
        a6.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
        a6.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)

        
        fig.set_facecolor("black")
        # plt.title(str(zz))
        if zz < 10: 
            strzz='0'+str(zz)
        else:
            strzz=str(zz)
        
        
        plt.savefig(savedir+'saturn_movie'+strzz+'.pdf', dpi=300)
        plt.savefig(savedir+'png_saturn_movie'+strzz+'.png', dpi=300)

        
        plt.show()

        naxis_mapcount[:]=0
        naxis_h3p[:]=0
        naxis_methane[:]=0
        naxis_heat[:]=0
        naxis_reflect[:]=0
        naxis_F323[:]=0
        naxis_ra[:]=0
        naxis_dec[:]=0
         
        naxis_mapcount[yc,xc]=m_counts
        naxis_h3p[yc,xc]=m_h3p
        naxis_methane[yc,xc]=m_methane
        naxis_heat[yc,xc]=m_heat
        naxis_reflect[yc,xc]=m_reflect
        naxis_F323[yc,xc]=m_F323
        naxis_ra[yc,xc]=m_ra
        naxis_dec[yc,xc]=m_dec
 
        start_time = geo.obs_start
        cml = (geo.lon_sun+360)%360
        geofirst=geo

        
    # if plot:

    # _plot(px, py, xc, yc, area, slices)
zz=zz+1

# image_g = np.nan_to_num(naxis_h3p/naxis_mapcount)
# image_g=np.nan_to_num(image_g)**0.5
# # print(np.nanmax(image_g))
# # image_g=image_g/np.nanmax(image_g)
# image_g=image_g/135
# image_r = np.nan_to_num(naxis_heat/naxis_mapcount)  # /np.max(im_methane)
# # print(np.nanmax(image_r))
# # image_r=image_r/np.max(image_r)
# image_r=image_r/140000
# image_b = np.nan_to_num(naxis_F323/naxis_mapcount)  # /np.max(im_reflect)
# # print(np.nanmax(image_b))
# # image_b=image_b/np.max(image_b)
# image_b=image_b/90000
# image = make_lupton_rgb(image_r*5, image_g*5, image_b*10)
        

# for xxx in range(180): aaa[:,xxx]=aaa[:,xxx]/np.nanmax(aaa[:,xxx])



image_g = np.nan_to_num(naxis_h3p/naxis_mapcount)
image_g=np.nan_to_num(image_g)**0.5
image_g=image_g/135*3
image_g[image_g != 0]= image_g[image_g != 0]-0.2
image_g[image_g < 0]=0

image_r = np.nan_to_num(naxis_heat/naxis_mapcount)  # /np.max(im_methane)
# print(np.nanmax(image_r))
# image_r=image_r/np.max(image_r)
# image_r=image_r/140000*1.1
image_r=image_r/140000*3
image_F323 = np.nan_to_num(naxis_F323/naxis_mapcount)  # /np.max(im_reflect)
image_reflect = np.nan_to_num(naxis_reflect/naxis_mapcount)  # /np.max(im_reflect)
image_b = np.nan_to_num(naxis_methane/naxis_mapcount)  # /np.max(im_reflect)
# print(np.nanmax(image_b))
# image_b=image_b/np.max(image_b)
# image_b=image_b/90000*5
image_b=image_b/90000*150
image = make_lupton_rgb(image_r*5, image_g*5, image_b*10)


print(zz,np.nanmax(naxis_lon))

# naxis_lon=naxis_lon2
# image_lat=np.nan_to_num(naxis_lat)
# image_lon=np.nan_to_num(naxis_lon)

image_ra=np.nan_to_num(naxis_ra/naxis_mapcount)
image_dec=np.nan_to_num(naxis_dec/naxis_mapcount)

image_ra2=  ((image_ra-(x_offset*xscale))/xscale/-3600.0)+geofirst.ra_target
image_dec2= ((image_dec-(y_offset*yscale) )/yscale/3600.0)+geofirst.dec_target

image_ra2[image_ra == 0]=0
image_dec2[image_dec == 0]=0

        

image_lat=np.zeros_like(image_ra)*np.nan
image_lon=np.zeros_like(image_ra)*np.nan


for xx in range(image_ra[:,0].size):
    for yy in range(image_ra[0,:].size):
    
        
        if image_ra2[xx,yy] != 0:
            geopixel=geofirst.pixel_params(image_ra2[xx,yy],image_dec2[xx,yy])
            image_lat[xx,yy]=geopixel['lat']
            image_lon[xx,yy]=geopixel['lon']
   
# if image_ra2[xx,yy] != 0:
#     rara=((yy-(x_offset*xscale))/xscale/-3600.0)+geofirst.ra_target
#     decdec=((xx-(y_offset*yscale) )/yscale/3600.0)+geofirst.dec_target
#     geopixel=geofirst.pixel_params(rara,decdec)
#     image_lat[xx,yy]=geopixel['lat']
#     image_lon[xx,yy]=geopixel['lon']

if np.abs(geofirst.lon_sub) < 90:
    image_lon2=image_lon+0.
    image_lon2[image_lon2>180]=image_lon2[image_lon2>180]-360
    image_lon2[np.nonzero(image_lon2)]=image_lon2[np.nonzero(image_lon2)]+180
    image_lon=image_lon2+0.


# image_lat2=image_lat+0
# image_lon2=image_lon+0
# for xxx in range(180): aaa[:,xxx]=aaa[:,xxx]/np.nanmax(aaa[:,xxx])

        
fig = plt.figure(figsize=(16,9))

scaling = ((700*3/16)/(460*2/9))
# xwidth*xnum*16/yheight*ynum*9

mytext=start_time[0:10]+' '+start_time[11:19]

a0 = plt.subplot(position=[0,1.0/scaling,1,0.164])
a0.set_axis_off()
a0.text(0.5, 0.7, mytext,c='w',stretch='ultra-expanded',family='Arial Black', ha="center", va="center", weight='bold',size=24)

a1 = plt.subplot(position=[0,0.5/scaling,1/3,0.5])
# a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
a1.set_axis_off()        

a1.imshow(image,aspect=1,origin='lower')
# plt.title(str(i))
# plt.plot([0,180],[550,550],linestyle='--',color='grey')
a1.set_xlim((1*xscale,8*xscale))
a1.set_ylim((0.9*yscale,5.5*yscale))

a1.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
a1.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
a1.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)

a2 = plt.subplot(position=[1/3,0.5/scaling,1/3,0.5])
# a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
a2.set_axis_off()        

a2.imshow(image_g,aspect=1,origin='lower',vmin=0,cmap=mygreen)
# plt.title(str(i))
# plt.plot([0,180],[550,550],linestyle='--',color='grey')
a2.set_xlim((1*xscale,8*xscale))
a2.set_ylim((0.9*yscale,5.5*yscale))

a2.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
a2.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
a2.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)


a3 = plt.subplot(position=[2/3,0.5/scaling,1/3,0.5])
# a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
a3.set_axis_off()        

a3.imshow(image_r,aspect=1,origin='lower',vmin=0,cmap=myred)
# plt.title(str(i))
# plt.plot([0,180],[550,550],linestyle='--',color='grey')
a3.set_xlim((1*xscale,8*xscale))
a3.set_ylim((0.9*yscale,5.5*yscale))

a3.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
a3.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
a3.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)


a4 = plt.subplot(position=[0,0.0,1/3,0.5])
# a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
a4.set_axis_off()        

a4.imshow(image_reflect,aspect=1,origin='lower',vmin=0,cmap='Greys_r')
# plt.title(str(i))
# plt.plot([0,180],[550,550],linestyle='--',color='grey')
a4.set_xlim((1*xscale,8*xscale))
a4.set_ylim((0.9*yscale,5.5*yscale))

a4.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
a4.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
a4.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)

a5 = plt.subplot(position=[1/3,0.0,1/3,0.5])
# a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
a5.set_axis_off()        

a5.imshow(image_b,aspect=1,origin='lower',vmin=0,cmap=myblue)
# plt.title(str(i))
# plt.plot([0,180],[550,550],linestyle='--',color='grey')
a5.set_xlim((1*xscale,8*xscale))
a5.set_ylim((0.9*yscale,5.5*yscale))


a5.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
a5.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
a5.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)

a6 = plt.subplot(position=[2/3,0.0,1/3,0.5])
# a31i=a3x1.contourf(lon2d, lat2d, mapmap_scaled2,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot',vmax=maxint,vmin=minint) 
a6.set_axis_off()        

a6.imshow(image_F323,aspect=1,origin='lower',vmin=0,cmap='copper')
# plt.title(str(i))
# plt.plot([0,180],[550,550],linestyle='--',color='grey')
a6.set_xlim((1*xscale,8*xscale))
a6.set_ylim((0.9*yscale,5.5*yscale))

a6.contour(image_lon,levels=minorlontick,colors='w',linestyles='dotted',alpha=0.4)
a6.contour(image_lon,levels=majorlontick,colors='w',linestyles='dotted',alpha=0.9)
a6.contour(image_lat,levels=np.linspace(0,90,num=int(90/15)),colors='w',linestyles='dotted',alpha=0.4)


fig.set_facecolor("black")
# plt.title(str(zz))
if zz < 10: 
    strzz='0'+str(zz)
else:
    strzz=str(zz)


plt.savefig(savedir+'saturn_movie'+strzz+'.pdf', dpi=300)
plt.savefig(savedir+'png_saturn_movie'+strzz+'.png', dpi=300)


plt.show()

'''
ffmpeg -framerate 4.348 -pattern_type glob -i '*.png' \
   -c:v libx264 -pix_fmt yuv420p out_1s100m.mp4